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By exploiting structural differences between electromagnetic and hadronic showers in a multivari- 
ate analysis we present an efficient Electron-Hadron discrimination algorithm for liquid argon time 
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I. INTRODUCTION 

Liquid Argon Time Projection Chambers (LAr-TPCs) 
are recognised as being a potential detector technol- 
ogy for a next-generation Neutrino oscillation experi- 
ment and their development is the subject of vibrant 
R&D programmes in Europe, Japan and the USA (see 
Ref. m and references therein). The simultaneous track- 
ing and calorimetry capability, with millimetric granu- 
larity, makes LAr-TPCs ideal to image Neutrino inter- 
actions over a wide range of energies as demonstrated 
recently by ICARUS 0. 

Despite the physics promise of LAr-TPCs, a fully au- 
tomated (that is, with zero human interaction) software 
system for reconstructing events has proven difficult to 
achieve. With the availability of more and more event 
data, recent years have seen progress in limited cluster- 
ing of structures and tracking. However, a general event 
is likely to contain a mixture of showers and tracks origi- 
nating from an unknown vertex anywhere in the sensitive 
volume. One major challenge for automated reconstruc- 
tion is to disentangle the fundamentally different shower 
and track structures and hence segment the event into its 
constituent final states. Accurate segmentation is critical 
for subsequent particle identification on the substructures 
and final classification of the interaction. Neutrino oscil- 
lation analyses proposed for next-generation facilities rely 
heavily on accurate interaction classification, motivating 
studies of this task. 

In this article, we concentrate on the particle identifi- 
cation step assuming exact segmentation, deferring dis- 
cussion of techniques for this more challenging task to 
later publications. We present a study using a set of 
effective and computationally simple variables which al- 
low the discrimination of Electron shower structures from 
all relevant hadronic shower structures in a fine-grained 
tracking detector such as a LAr-TPC. For the intended 
physics application, the dominant Hadron shower struc- 
tures result from Pious, Protons or Kaons in the final 
state. In this context, Electron-Hadron shower discrim- 
ination is important for Neutrino oscillation studies in 
order to measure Electron-neutrino events in a Muon- 
neutrino beam, resulting in an energetic Electron in the 
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LAr medium after a charged-current interaction with a 
nucleus. Therefore, the Electron is considered to be the 
"signal" and anything non-electron, i.e. Hadron showers, 
is labelled as "background" , unless stated otherwise. 

Additionally, the energy resolution of any LAr-TPC 
will depend on the nature of the particle depositing en- 
ergy due to quenching factors. Electromagnetic and 
hadronic particle energy depositions will each require 
separate energy calibrations. Stronger ionising parti- 
cles such as Hadrons would show particle specific rela- 
tive quenching factors (relative to Electrons), hence for 
their energy measurement particle identification would 
be beneficial. 

The set of variables also prove to be useful for other 
particle identification tasks. Track structures associated 
with Protons or Muons can be identified by these vari- 
ables as well as, to a lesser extent, hadronic final states 
such as neutral and charged Pious and Kaons. 

In the following, we describe the analysis algorithm by 
first introducing the discrimination variables and then 
presenting their applicability to identifying particles in 
various Neutrino interaction event topologies. 



II. ELECTRON-SHOWER DISCRIMINATION 
ALGORITHM 

Below, we detail the calculation and physics motiva- 
tion of six discriminating variables. The main assump- 
tion underpinning all of the following is that clustered 
energy deposits, hits, in the detector from a particular 
particle have been collected as a single object by the re- 
construction, i.e. a cluster. 

This crucial first step of any LAr detector analysis is as- 
sumed to have taken place already. It should be stressed 
that a reliable, complete and automatic initial clustering 
step is, in fact, among the unsolved challenges of LAr 
data analysis, owing to the richness of topological event 
structures in such a detector, and will not be addressed 
here. 

The first step of the discrimination algorithm consists 
of transforming the isolated cluster using a principal com- 
ponent analysis (PC A) [3]. The data at this stage con- 
sists purely of detector hits forming a point-cloud of un- 
known shape. Dealing with LAr detectors translates into 
the capability of assigning three spatial coordinates and 
a charge value to each individual hit (with corresponding 
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FIG. 1: The set of six discrimination variables applied to the CN2PY beam (see text) Neutrino events. Signal events are defined 
as electromagnetic shower events from Electrons, compared to the complete shower background, which includes all simulated 
structures from Protons, charged Pions, neutral Pious, charged and neutral Kaons and Muons, all with appropriate weighting 
factors. The vertical text on the right-hand side indicates the under- and overflows 0|. 
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FIG. 2; Similar to Fig. 1; the set of six discrimination variables here applied to uniform energy final state particles, sampled 
from 10 MeV to 4.5 GeV. 



uncertainties) . 

The reason for the initial PCA step is that it allows us 
to consider the unknown point-cloud to be aligned and or- 
dered in a consistent way. The PCA is set up such that it 
calculates the three principal spatial components (ignor- 
ing the charge of the hit) as eigenvalues of the transfor- 
mation matrix, which itself consists of the corresponding 
eigenvectors. Applying the transformation aligns the co- 
ordinate axes to the principal axes. This allows to utilise 



transverse and longitudinal geometric information due to 
the presence of a well defined axis for the structure. 



A. The lateral projection variable "lat" 

This discrimination variable is motivated by previous 
and ongoing accelerator experiments and their analysis of 
electromagnetic (EM) and Pion showers (see for example 
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FIG. 3: Similar to Fig. 1; the set of six discrimination variables here applied to uniform energy final state particles, sampled 
from 10 MeV to 770 MeV. 



Rcf. "s"]). The entire point-cloud is projected onto the 
plane of minor principal components (the lateral plane 
to the major principal axis). A core region is defined 
around the projected mean value and the ratio of energy 
(charge) deposited in the core to the total energy (charge) 
yields the "lat" variable. 

Charged Hadron showers tend to have a denser core 
region compared to EM showers, or often resemble a 
track- like structure with only a thin "halo" around it, or 
dissolve into a wide-spread cloud of hits. Proton tracks 
might or might not shower at all. EM showers, by con- 
trast, tend to produce more constant halo-to-core ratios, 
as shown in Fig. [TJ 

The precise definition of the core region turns out to 
be non-critical. For simplicity, a circular shape is chosen, 
which is centred on the mean value of the projected point- 
cloud and is defined with the Moliere radius in liquid 
argon, which is equal to 9.61cm Q. Each hit inside the 
core contributes to the sum of energy in the core, which 
is then divided by the total energy of the shower (the 
sum of the energies of all hits). 



minimum threshold radius is set in order to avoid the 
Coulomb singularity (for practical purposes chosen as 
1 mm in this analysis) , and finally the hit concentration 
is calculated as the sum of charge (energy) of the hits 
divided by their radial distance. 

The idea is to calculate the Coulomb potential energy 
of the point-cloud around the main principal axis. This 
approach delivers complementary information to the pre- 
vious rigid division into two parts, since it emphasises the 
extremes such as track-like structures as being particu- 
larly highly concentrated, and very loose structures as be- 
ing particularly low in concentration. Figures [U [2] and [3] 
show the main benefit of this variable. The Hadron show- 
ers are dominated by less densely concentrated structures 
(they peak at low Coulomb energies around the main 
principal axis). A second population of Hadron events 
contain Proton tracks as well as sufficiently dense Pion 
tracks, contributing a large value to this concentration 
variable. EM showers are rather non-specific but are sig- 
nificantly more concentrated than Hadron showers. 



B. The hit concentration variable "con" 

While the previous variable emphasises a subdivision 
of the shower structure into two main components (core 
and halo), the hit concentration variable "con" attempts 
to use the same lateral information globally. This is in- 
dependent of any material constant, such as the Moliere 
radius. 

The variable first calculates the radial distance of each 
hit from the main principal axis, i.e. the radial dis- 
tance from the origin in the lateral projection plane. A 



C. Initial ionisation loss variable "dedx" 

The above two variables already provide excellent 
EM/Hadron shower discrimination potential. This third 
variable is mainly introduced to improve the low energy 
discrimination of the analysis, enabling additional appli- 
cations of the procedure. 

This variable attempts to use the information about 
the longitudinal shower characteristics. At Hadron col- 
liders the entire longitudinal shower profile is used @, 
but it was found that only the first part of any shower 
is necessary for providing discrimination at low energies. 
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This variable is ineffective at high primary particle ener- 
gies and could be dropped since the ionisation loss curves 
of final state particles become indistinguishable. 

The calculation proceeds as follows: partition the en- 
tire set of hits into roughly equal parts in hit numbers and 
consider only the first slice. This slicing is particularly 
simple after the PCA and hit ordering. The first slice 
then contains hits from the start of the shower. Then, 
calculate the sum of charges (energies) contained in the 
slice, its physical length along the major principal axis 
and build the specific energy loss variable dE/dx, nor- 
malised to the number of hits in the slice. For quantita- 
tive comparisons, energies were normalised to the liquid 
argon critical energy of 30.6 MeV and lengths to the spe- 
cific radiation length of 14 cm [gJ. 

The number of slices is fixed to 50, since fewer slices 
provide less discriminating power, while too many slices 
diminish the number of available hits in each slice, en- 
hancing statistical fluctuations. Typically, each slice 
would contain more than 40 hits, sometimes increasing 
to several hundred hits for large showers, depending on 
the particle energy. 




Energy [GeV] 



FIG. 4: The normalised CN2PY Neutrino beam fluxes ^ 
from negative Pion decays. The black line shows the Muon- 
neutrino flux; red - the Electron-neutrino flux; green - anti- 
Muon-neutrino flux and blue - anti-Electron-neutrino flux. 

As expected, this variable delivers additional discrimi- 
nation power against low energy Protons with relatively 
large specific ionisation energy loss compared to Pious, 
Electrons and Muons. It might also support distinguish- 
ing low energy heavy charged fragments and their cor- 
responding showers. At high primary Neutrino energies, 
this variable becomes ineffective, since the resulting Pro- 
ton energies are generally high enough not to be easily 
distinguishable according to their ionisation energy loss 
(and Pious even less so). Note that the discrimination 
power of this variable is not clearly visible in Fig. [31 Due 
to the low energies of final state particles in this case, 
the calculation would fail (too few hits in a slice) and 
the variable does not contribute to the analysis in these 



cases. The rigid definition from above would have to be 
optimised further for such a low energy study. 



D. Spatial extent variables "extx,y,z" 

The final three proposed discrimination variables are 
obtained by quantifying the spatial extent of the hit cloud 
in the three available principal component axes. A def- 
inition of the extent in purely geometrical terms turns 
out to be very useful. 

The spatial extent is calculated by obtaining the con- 
vex hull in 3D around the hit cloud and finding the 
maximum distance (extent) of the hull on all three axes. 
Building the convex hull around an object is a standard 
algorithm in image analysis . Even though its applica- 
tion in 3D is rather exotic, ready-made implementations 
of the code exist, for instance in the scientific python li- 
brary, SciPy fs]. The meaning of the term convex hull 
can be visualised by considering stiff packaging paper, 
wrapped around an irregularly shaped object, i.e. the 
hit cloud, in order to enclose it completely. The paper 
would represent the convex hull, without the folds and 
overlaps real paper would produce. 

The convex hull is obtained by calculating the De- 
launay tessellation Q of the hit cloud and finding the 
outermost Delaunay triangles, made of data members, 
enveloping the entire cloud. The difference between the 
maximum and minimum distance on each (principal com- 
ponent) coordinate axes is calculated, giving the extent 
of the hull in all three dimensions. 



III. RESULTS 

Several simulated single particle event classes were pre- 
pared as input to the analysis. We have tested the al- 
gorithm with a Muon-neutrino flux generated as part 
of the LAGUNA-LBNO FP7 programme ^ to evaluate 
the physics potential of a long baseline experiment be- 
tween CERN and an underground far detector location 
in Pyhasalmi Finland (CN2PY). The CN2PY Neutrino 
flux jlO| (un-oscillated) from decaying negatively charged 
Pions (for a distance of 2300 km) provided the input for 
the first case study. Neutrino reaction final states were 
produced using the GENIE event generator JTj. The 
normalised Neutrino energy distributions are shown in 
Fig. m with relative weighting factors given in Table H] 

TABLE I: The relative weighting factors for the CN2PY 
Neutrino beam (un-oscillated) from negatively charged Pions 
taken as input for this study. 

I'fi flux flux i/e flux l/e fluX 

1.0 3.471 X 10'^ 4.455 x 10~^ 1.0772 x 10"^ 
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The GENIE simulated events were un-filtered, with ah 
possible final state particles allowed. The total number of 
each individual particle final state from GENIE in each of 
the four Neutrino fluxes then results in individual weight- 
ing factors for that particle species during the subsequent 
analysis. The complete set of weighting factors is listed 
in Table |TT1 Final state Hadrons which do not decay by 
default in Geant4 have been ignored in the subsequent 
analysis. 

These single particle energy distributions from the GE- 
NIE simulations then serve as the particle source in the 
Geant4 Monte-Carlo transport simulation toolkit p^ . 
The simulation models an homogeneous liquid argon de- 
tector as a cylinder with a radius and height of 100 m, 
in order to fully contain the showers. Particles start 
at the centre of the cylindrical volume and are tracked 
through the volume with interactions modelled by the 
QGSP J3IC_IIP physics list, with all electromagnetic and 
hadronic processes enabled. The LAr volume is divided 
into (1x1x1) mm"^ cubic voxels, and all primary and 
secondary particles are tracked through voxels down to 
zero energy, or until they leave the volume. Energy de- 
posits by charged particles passing through voxels are 
tallied into a map between voxel coordinates and the to- 
tal energy deposited. No attempt was made to model the 
readout system, but quenching of ionisation events spe- 
cific to LAr detectors was taken into account following 

El. 

Results are collected and processed in a final step 
by four selected multivariate analysis algorithms (MVA) 
taken from the TMVA toolkit Each MVA is based on 
supervised learning, requiring separate training and test 
data samples. Each sample is chosen from a random split 
in half of the available data. The final evaluation delivers 
optimal cuts for each selected MVA method, producing 
efficiencies for signal selection, background selection and 
corresponding purities. 

Our second case study attempts a pure algorithm test 
by running the identical analysis as discussed above over 
a sample of final state particles sampled from individ- 
ual uniform energy distributions. The high energy case 
(labelled Uniform HE in Tab. IIIip samples particle en- 
ergies between lOMeV (typical liquid argon threshold 
energy value) and 4500 MeV (approximately the mean 
Neutrino energy of the proposed CN2PY experiment), 
all with equal weighting factors, while the low energy 
option (labelled Uniform LE in Tab. IIIII) samples up to 
770 MeV (corresponding to the mean energy of the cur- 
rent JPARC neutrino beam). These case studies gauge 
the sensitivity of the method to knowledge of the beam 
kinematics. The proposed variables are independent of 
beam characteristics but the multivariate analysis step 
depends on these case study assumptions. 



A. Electron to combined Hadron/lepton 
background discrimination 

Results in this section summarise the performance of 
the proposed discriminating variables when attempting 
to identify Electrons in a data sample consisting of all 
dominant final state particles, i.e. charged and neutral 
Pion, Protons, Kaons and Muons. The most straightfor- 
ward way of illustrating the effectiveness of the proposed 
variables can be seen in Fig. [5l showing the background 
rejection as a function of signal efficiency for all four cho- 
sen methods, the nearest neighbour (KNN), two varia- 
tions of boosted decision trees (BDT, BDTG) and a neu- 
ral network (MLPBNN). These plots would ideally show 
a complete background rejection for all signal efficien- 
cies up to 100%. The closer the curves follow these ideal 
lines on the figure, the better the discrimination of the 
selected MVA method. The MVA methods mentioned 
above have been chosen as the most effective among a 
list of 24 possible methods offered in the TMVA toolkit. 

Quantitative results for efficiencies and purities in EM 
shower discrimination are given in Table IIIII and illus- 
trate a successful (well above 90% efficiency and purity) 
analysis technique in distinguishing EM showers from 
Hadron showers in a LAr detector at all interesting Neu- 
trino energies. Additionally, we note that the reduced 
efficiency and purity values for low energy primary neu- 
trinos compared to the high energy neutrino case study 
is as expected. Finally, the method appears relatively 
insensitive to presence or absence of priori beam knowl- 
edge. 



B. Additional applications 

The variables discussed above lend themselves to fur- 
ther applications simply because they are sensitive to 
only the geometrical shapes of clouds of hits in 3D space 
as opposed to kinematic variables. As a first application 
we consider neutral current (NC) interactions in LAr. 
NC interactions form an important background to Neu- 
trino oscillation signals. The dominant detectable sig- 
nature of such interactions would be the appearance of 
neutral Pions in an event. Due to their prompt decay into 
a pair of photons, with the inevitable EM showers, they 
constitute the most destructive background to Electron- 
neutrino charged-current interactions. This means that 
identifying final state tt*^ would be beneficial for several 
reconstruction tasks. One study [Ti^ delivered a strong 
discrimination factor of 10^^ for mis-identified tt*^ in a 
LArTPC, a number which is widely quoted. Our efforts 
here deliver a weaker 7r° discrimination, see below. How- 
ever it does not rely on complex computational tasks such 
as shower separation and fitting and subsequent vertex 
finding in a potentially very busy event structure, i.e. the 
reason why 14] relies on visual event selection. 

The optimal signal and background efficiencies for the 
7r° selection are shown in Table IIIII (see also Fig. [5] and 
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TABLE II: The relative weighting factors for the final state particles resulting from the CN2PY Neutrino beam containing 
particles given in the first column. These are each taken into account during the multivariate analysis stage. Note that Muon 
and Electron energy spectra are taken from Genie simulations using the corresponding neutrino-flavour weighting factor. 





7T+ 




7V° 


P 


A-± 




7 




0.203 


0.508 


0.409 


1.610 


0.016 


0.022 


0.035 




2.42 X 10"^ 


9.89 X 10"^ 


2.04 X 10"^ 


6.284 X 10"^ 


1.43 X 10"^ 


1.1 X 10"^ 


1.96 X 10"^ 




9.4 X 10"* 


2.32 X 10"^ 


1.85 X 10"^ 


7.2 X 10"^ 


7.7 X 10"^ 


1.1 X 10~* 


1.8 X 10"* 




7.6 X 10"* 


3.1 X 10~* 


6.4 X 10~* 


1.95 X lO'^* 


4.5 X 10"^ 


3.4 X 10"^ 


5.7 X lO"'' 



g 0.£ 

■a 0.£ 
c 

S 0.7 

Ui 

0.£ 

0.; 
0.: 



Ill 


III.. 


~ • • ] 
























































M\ 




od: 


















BDT 




















BDTG 


















KNN 








, , , 









0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Signal efficiency 



Background rejection versus Signal efficiency 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 



JIM 














1 




































































MV/a MothnH* 




















BDT 
















E 






KNN 




















BDTG 

















0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Signal efficiency 



FIG. 5: Background rejection as a function of the signal efficiency for all four chosen methods for a (left panel) high energy 
uniform Electron energy distribution and (right panel) CN2PY beam. The methods are: the k-nearest neighbour search 
(KNN), two variations of boosted decision trees (BDT, BDTG) and a neural network (MLPBNN), all originating from the 
TMVA toolkit _4]. For numerical results on this analysis, see Table Hill 
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FIG. 6: Variables for a tt final state particle signal from a CN2PY beam compared to the complete shower background. The 
final three variables describing the geometrical shape of the structures are improving the discrimination when compared to a 
three (the first three) variable analysis. 



Fig. [7|)- Just as for the study on Electron to Hadron discrimination, all subsequent studies on changing the 
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FIG. 7: Variables for a vr'' final state particle signal from a uniform energy distribution compared to the complete shower 
background. 



signal, here to the tt'^ particle, consider the complete set 
of background particles, listed as final state particles in 
Table im and additionally the relevant leptons, i.e. Muons 
and Electrons. All final state particles other than the 
selected signal particle constitute background. 

Other final state particles can also be identified with 
varying degrees of success. In Table HIl we list results for 
Muon, Protons, Kaons, where charged Pion results are 
practically identical to Kaon results, hence have been left 
out. 



C. Detector threshold effects 

Detector-specific effects depend heavily on our primary 
assumption of successful clustering of hits belonging to 
a single particle in an overall unspecified detector event. 
However, any automated analysis will at some point have 
to solve this challenge for finely-grained LAr detector de- 
vices. 

Nevertheless, threshold variations would remove (or 
add) hits to the total simulated single particle event and 
potentially change geometrical structures. In order to 
test the robustness of our method against such a real- 
istic detector effect, we studied two scenarios: one real- 
istic hit energy cut at 120keV/mm for quenched (mea- 
sured) energies and one essentially without an hit energy 
cut (IkeV/mm). The threshold choice of 120keV/mm 
originates from the simulated onset of the Muon energy 
loss peak which would be left practically fully intact by 
this cut. The minimal cut scenario at 1 keV/mm roughly 
doubles the total amount of hits for most single particles 
simulated. 

The change from these threshold effects on all of the 



particle identification results shown earlier is negligible, 
typically less than 0.3% on any signal or background ef- 
ficiency quantity. 

Similar robustness should be pointed out with respect 
to hit smearing effects, i.e. the case when individual hits 
are represented by more realistic volumes of three to six 
mm cubes (or even asymmetric dimensions). Our pro- 
posed variables are purely measures of the geometry of 
structures in a LArTPC. As long as the smallest spatial 
measure, i.e. a hit, is small compared with the extent of 
shower/track structures, we do do not expect any of our 
analysis results to change significantly. 



IV. CONCLUSION 

We have presented a study identifying particles in 
event structures obtained in high-granularity detectors 
such as a liquid argon TPC using six discriminating vari- 
ables. Combining them using MVA algorithms from the 
TMVA toolkit gives signal efficiencies for Electron and 
Muon identification well above 90% with high purities 
on all Neutrino interaction channels. Protons, Kaons and 
charged Pious (not discussed above due to very similar 
results to Kaons) show anything from acceptable to poor 
purity but high efficiency. Neutral Pion identification 
can be made very efficient if a purity between 16% and 
25%, depending on the initial Neutrino energy, is not a 
concern. 

The data used in this work has been assumed to con- 
sist of hits containing three spatial coordinate values and 
one energy value, which should be the case for LAr TPC 
detectors. Furthermore, it is assumed that each single 
structure to be analysed has been clustered correctly 
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TABLE III: Particle identification efficiencies and contamination (background efficiency or 'contamination' is equal to 1.0 - 
signal purity) in a combined background of Hadron showers and Muons for all selected MVA methods. Efficiencies are calculated 
for the optimal decision boundary values defined by maximising the signal-to-total event ratio, S/ {S + B). 
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FIG. 8: Similar to Fig. 1 for the CN2PY beam producing Muons in the final state, here defined as signal compared to the 
complete shower background. 



which is still an outstanding challenge for automated 
data analysis on high granularity detector data. Nev- 
ertheless, if such a clustering can be achieved by some 
means, then this work proposes an efhcient and reli- 
able method to identify Electron (as well as positron) 
events in a combined hadronic and leptonic background. 
These would be the primary signatures of charged cur- 



rent Electron-Neutrino interactions which constitute the 
main target of future long baseline experiments. It also 
allows the identification of neutral Pions without addi- 
tional knowledge of decay vertices (or similar quantities) 
in a mixed background of other Hadrons and leptons, and 
can be used for Muon identification and as an efficient 
first pass on Proton, Kaon and charged Pion identifica- 
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FIG. 9: Similar to Fig. 1 for the CN2PY beam producing Protons in the final state, here defined as signal compared to the 
complete shower background. 
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FIG. 10: Similar to Fig. 1 for the CN2PY beam producing Kaons (charged and neutral combined) in the final state, here 
defined as signal compared to the complete shower background. 
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